Dynamic Pathways for Viral Capsid Assembly 
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We develop a class of models with which we simulate the assembly of particles into Tl capsid-like objects 
using Newtonian dynamics. By simulating assembly for many different values of system parameters, we vary 
the forces that drive assembly. For some ranges of parameters, assembly is facile, while for others, assembly 
is dynamically frustrated by kinetic traps corresponding to malformed or incompletely formed capsids. Our 
simulations sample many independent trajectories at various capsomer concentrations, allowing for statistically 
meaningful conclusions. Depending on subunit (i.e., capsomer) geometries, successful assembly proceeds by 
several mechanisms involving binding of intermediates of various sizes. We discuss the relationship between 
these mechanisms and experimental evaluations of capsid assembly processes. 



I. INTRODUCTION 

This paper is devoted to introducing a simple class of cap- 
somer models, and demonstrating that Newtonian dynamics 
of these models exhibit spontaneous assembly into 60-unit 
icosahedral capsids, depending upon conditions (i.e., particle 
concentration and force field parameters). We believe it is the 
first report of statistically meaningful simulations of capsid 
assembly that follow from unbiased dynamics obeying time- 
reversal symmetry and detailed balance. 

The formation of viral capsids is a marvel of natural engi- 
neering and design. A large number (from 60 to thousands) of 
protein subunits assemble into complete, reproducible struc- 
tures under a variety of conditions while avoiding kinetic and 
thermodynamic traps. Understanding the features of capsid 
components that enable such robust assembly could be impor- 
tant for the development of synthetic supra-nano assemblies. 
In addition, this knowledge is essential for the development 
of anti-viral drugs that inhibit capsid assembly or disassembly 
and could focus efforts to direct the making of highly specific 
drug delivery vehicles. These goals necessitate the ability to 
manipulate when and where capsids assemble and disassem- 
ble. Thus, we seek to determine what externally or internally 
controlled factors promote or alleviate dynamic frustration in 
the capsid assembly process. Although many viruses assem- 
ble with the aid of nucleic acids and scaffolding proteins, the 
first step toward this objective is to understand the inherent 
ability of subunit-subunit interactions to direct spontaneous 
assembly. 

The equilibrium properties of viral capsids have been the 
subject of insightful theoretical investigations (e.g Refs. QIS 
|^|j,|5l|^ and the assembly process has been investigated in a 
number of experiments (e.g. Refs. 7, 8, 9, 10, 11, 12.mfl4l 
O, yet this process is still poorly understood for many viruses 
(e.g. Ref. not) . Assembly is difficult to analyze experimen- 
tally because most intermediates are transient. With single 
molecule techniques, it is now possible to directly probe inter- 
mediate structures. Each intermediate, however, is a member 
of a large ensemble of structures and pathways that comprise 
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the overall assembly process. Formation of an intermediate re- 
quires collective binding events that are regulated by a tightly 
balanced competition of forces between individual subunits. It 
is difficult, with experiments alone, to parse these interactions 
for the factors that are critical to the assembly process. Thus, 
it is useful to have complementary computational models in 
which the effects of different interactions can be isolated and 
monitored. 

Studying assembly through computation is challenging be- 
cause short range subunit-subunit properties regulate the for- 
mation of overall structure. Binding and unbinding rates of 
individual subunits are orders of magnitude faster than the 
overall assembly times. Furthermore, these rates are con- 
trolled by interactions defined on atomic lengths, which are 
three orders of magnitude smaller than typical capsid sizes. 
Prior computational studies provide valuable insights that we 
build upon; in particular, Zlotnick pioneered a rate-equation 
description of assembly | ij] and Berger and co-workers de- 
veloped particle based methods |fl8|. Earlier studies, although 
an important foundation for our work, are limited in that they 
have been based upon pre-conceived pathways of assembly 
LI 7.. ,19.. ,20., .2L i22>]. or dynamics that did not obey detailed 
balance Il8l l23l I25I1 . or dynamics that was anecdotal 
26, 27, 28"]. These approaches can be useful and physi- 
cal justifications for them can be made. Nevertheless, we seek 
to avoid these limitations in order to understand the nature of 
possible kinetic traps and the extent of ensembles of success- 
ful assembly events. In Section |II| we present our class of 
models for capsid subunits. We evaluate the thermodynamic 
properties of this model in Section |III| and then discuss the 
results of dynamical simulations in Section IIVI By simulat- 
ing assembly for many different values of system parameters, 
we vary the strength of the forces that either drive or thwart 
assembly. We identify regions of parameter space in which 
two forms of kinetic traps prevail and we elucidate processes 
by which dynamic frustration is avoided in other regions of 
parameter space. 



II. MODEL 

Capsomers. Capsid proteins typically have several hun- 
dred residues that fold into well defined shapes with specific 
interactions that lead to attractions between complementary 
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sides of nearby subunits. We imagine that by integrating over 
degrees of freedom, such as atomic coordinates, as capsid 
proteins fluctuate about their native states, one can arrive at 
a model in which subunits have excluded volume and asym- 
metric pairwise bonding interactions between complementary 
sides. Several models have been presented in which asym- 
metric subunits, the capsomers, are represented by conglom- 
erates of spherically symmetric particles with varying interac- 
tion strengths |6, 23l. E7lE3l . These approaches can describe 
complex excluded volume shapes. The approach we take here, 
however, is simpler, and motivated by the modeling described 
in Ref. 26. Specifically, we use only a single spherical ex- 
cluded volume per capsomer, and we use internal bond vectors 
to capture the effects of protein shape and complementarity. 

Our model capsomers are of three types: B3, B4 and B5, 
which contain three, four and five internal capsomer bond vec- 
tors, respectively. These bond vectors, b,-"', are pictured in 
Figure 1 . The index a goes from 1 to rib, where n\, = 3 for 
the B3 model, nb = 4 for the B4 model, and so forth, and 
the index i goes from 1 to iV, where iV is the number of cap- 
somers in the system. The vector r|"' = Rj + b-"' is the 
position of interaction site a on capsomer i, where is the 
center of capsomer i. All bond vectors have the same mag- 
nitude, h. Within a capsomer frame of reference, the bond 
vectors are fixed rigidly. They move only because the cap- 
somer translates and rotates. This is not to say that proteins 
do not fluctuate. Those fluctuations, we imagine, have been 
averaged over, i.e., integrated out of the model at the level we 
consider. 




FIG. 1: Geometry of bond vectors in the B3, B4, and B5 capsomer 
models, the center of capsomer i is at Ri. The angles between indi- 
cated bond vectors within a capsomer are specified in degrees. They 
do not sum to 360° = 2-k because the bond vectors are not coplanar. 



The net potential energy of interaction among iV cap- 



somers, 2, . . . . N), is taken to be pair-decomposable, 

N 

C/(l,2,...,iV)= ^ u(i,2). (1) 

where u{i,j) depends upon the bond vectors and centers of 
capsomers i and j. The particular form for this pair poten- 
tial depends upon which of the three models, B3, B4 or B5 
is under consideration. In each case, however, the potential 
is constructed so that the lowest energy configurations coin- 
cide with separate icosohedral clusters of 60 identical cap- 
somers. These clusters represent capsids with triangulation 
number (T) of one |2]; for example, design B5 is consistent 
with Fig. 3 of Ref 29. Our B3 model is similar to the model 
considered in Ref. 26. 

In each model, bond vectors or interaction sites have com- 
plementary counterparts. For example, in the B3 model, in- 
teraction site pairs (a, f3) = (1, 2) and (3,3) are the primary 
complementary pairs. This means that a favorable potential 
energy of interaction between a pair i and j of B3 capsomers 
has two ways of occuring: 1) interaction site 1 on one cap- 
somer overlaps with interaction site 2 on the other capsomer, 
and the respective bond vectors b,-^^ and b^^^ are nearly an- 
tiparrallel; 2) interaction site 3 on one capsomer overlaps with 
interaction site 3 on the other, and b^- and b^ are nearly 
antiparrallel. The only favorable (i.e., attractive) interactions 
are those associated with primary complementary pairs. 

In addition, there are secondary complementary pairs. For 
example, in the B3 model, with primary complementary pair 
(1,2), there is the secondary pair (3,3). This means that a fa- 
vorable interaction affected by the primary complementary 
pair (1,2), as described in the previous paragraph, also re- 

(3) (3) 

quires that b,^ ' and b^' are nearly coplanar Similarly, for 
the primary complementary pair (3,3), the secondary pair is 

(3) (3) 

either (1,2) or (2,1), meaning that if b^ and b!j are antipar- 

allel, favorable interactions result only if b^-^^ and b^^^ are 

nearly co-planar and hf ^ and bj^-* are nearly co-planar. Of 

course, because the capsomers are rigid bodies, b|^' and b^^' 

being co-planar implies b^^ and b^^^ are co-planar The pri- 
mary and secondary pairs for each of the models are listed in 
the entrees to Table I. Local bonding associated with these 
complementarities and resulting capsid structures are illus- 
trated in Figure 2. In creating these pictures, it is imagined 
that excluded volume interactions prohibit an interaction site 
from participating simultaneously in more than one favorable 
complementary interaction, as is the case for the models we 
describe. 

The dependence of subunit-subunit interactions on the ori- 
entation of primary and secondary pairs incorporates the fact 
that there is a driving force for subunits to align complemen- 
tary regions to maximize the contact between complementary 
residues. Capsid curvature in the minimum energy orientation 
arises from that fact that the angles between bond vectors on 
a given subunit do not sum to 2tt. 

Pair potential. The potential energy of interaction between 
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TABLE I: Primary and secondary complementary pairs and associ- 
ated angles for the three capsomer models. 



two capsomers, say 1 and 2, is taken to have a spherically 
symmetric repulsive part, uo(|R2 — R-i|), and an attractive 
part that depends upon both R2 — Ri and the bond vectors 
associated with the two capsomers, 

u(l,2) = uo(|R2-Ri|) 

+ui(R2-Ri,{b(")},{bi^)}). (2) 

For the reptilsion, we have chosen the Weeks-Chandler- 
Andersen 13011 potential, 

uo{R) = 4e - {a/Rf + 1/4] , R < 2^/^ 

= 0, R> 2^l^a. (3) 

For the attractions we have chosen 
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where the primed sum is over primary complementary pairs, 

Um{t) = 4£b 
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FIG. 2: Complementary pairs and bonding of capsomers. The first 
column specifies the model, the second illustrates the local bond- 
ing consistent with the complementary pairs of bond vectors, and 
the third illustrates the resulting complete capsid, with bonds depict- 
ing the attractive interactions resulting from complementary pairs. 
The pictures of complete capsids and all simulation snapshots shown 
in this work were generated in VMD |49]. The size of the spheres 
in these pictures has been reduced to aid visibility; parameters are 
chosen such that the minimum energy distance between neighboring 
capsomers is at the minimum in the WCA potential, Eq.|21 



which has its minimum value, —Eb, when the separation of 
complementary pair interaction sites is zero, and Sq/3(1, 2) is 
the switching function, given by 



Sa/3(1,2) 



1 



COS 



a/3,1) 
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a/3,2) 
2 



(6) 



for models B3 and B5, and by 

Sa/3(1,2) = 



1 r 

4 



cos I -kO 



g("/3) 



X 



COS 



(a/3,1) , , 
12 / fa 



1 



(7) 



for model B4. The angle variables used in these expressions 
are defined in Table I. Notice from that table, specifying a 
specific primary pair of complementary bonds prescribes spe- 
cific corresponding secondary pairs. The switching function 
goes smoothly from 1 to as the angle variables 0^''\ (j)^^''^'' 
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and (^^^'"^^ change from to 6'in , 0m and respectively. 
Increase of these maximum angles and increases the 
configuration space in which two nearby subunits attract each 
other, but also weakens the driving force toward the minimum 
energy orientation. 

Dynamical simulations Dynamical trajectories were cal- 
culated using Brownian dynamics, in which particle motions 
are calculated from Newton's laws with forces and torques 
arising from subunit-subunit interactions as well as drag and 
a random buffeting force due to the implicit solvent. We use 
the following coupled equations of motion 



8t, 



(8) 



where us is the angular velocity, the force is given by 

F, = (9) 
and the torque is given by 



(a) 



(a) 



(10) 



while (5Fi and are a random force and torque, with covari- 
ances given by 

((5F,(i)(5F,(i')) = U(t-t')5,,2k^Th 
{dT,{t)STj{t')) = I5{t - t')5,j2k^T (11) 

where 1 is the identity matrix. The friction coefficients for 
translation and rotation are 7 and 71 , respectively, and k^T is 
the thermal energy. 



Parameter 


Value 


Definition 




1 


WCA energy parameter, Eg.lsl 




9-22 


Attractive energy strength, Eg.lsl 


b/cr 


2-5/6 


Bond vector length 


7,-/70-2 


0.4 


Rotational friction coefficient, Eq.|8 


4>m (rad) 


3.14 


maximum dihedral angle, Eq.|S| 


(rad) 


0.1 - 3.0 


maximum bond angle, Eg.lSI 


L/a 


11 - 100 


Simulation boxsize 


N 


1000 


Number of subunits 


Co = Afo-^L-s 


0.001 - 0.75 


Concentration of subunits 


r^/a 


2.5 


Attractive energy cutoff distance 


St/to 


0.006 


Timestep 


tabs /to 


6 X 10^ 


Final observation time, 10* steps 



TABLE II: Parameter values used for dynamical simulations in this 
work, where a is the unit of length, ksT is the thermal energy, 7 is 
the translational friction constant, Eq.|8| and to = ja'^ / (ASksT) is 
the unit time. 

In our implementations of these equations, rigid body rota- 
tions were performed with quaternions 1 3 1 1 and rotational and 
translational displacements were calculated using the second 
order stochastic Runge-Kutta method 1 3Z i,33J . as described in 
Appendix A. Periodic boundary conditions were used to sim- 
ulate a bulk system. We employed reduced units in which the 
particle diameter a — 1, k^T is the unit of energy, and time 



is scaled by to = ^a'^ / (iSksT). Each trajectory considered 
N = 1000 subunits and ran for 10* steps, usually with a time 
step of 0.006. The values of all parameters used in this work 
are documented in Table II. If units of length, cr, and temper- 
ature, T, are chosen to be cr = 2 nm and T = 300 K, the final 
observation time after 10* steps is tobs — 227.5 /is, subunit 
concentrations, Cq range from 2.08 x lO^** — 0.156 mol/L, 
and binding energies, Sb, range from 5.4 — 13.2 kcal/mol. 




FIG. 3: Binding free energies for dimerization calculated from 
Eas. II4II5I nine) and Monte Carlo simulations (points). Free en- 
ergies are with reference to a standard state with volume fraction of 
1 and free rotational motion, and the maximum angle parameters, 
defined in Eq.|6| are 6m = 0m = 0.5. 
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FIG. 4: The thermodynamic critical subunit concentration for capsid 
formation, pcc, as calculated from Eas ll5lfT6l and ^| for design Ba 
and 0m ~ TT. Above these subunit concentrations, most subunits will 
be found in complete capsids at equilibrium. 



III. THERMODYNAMICS OF CAPSID ASSEMBLY 

The equilibrium concentrations of free subunits 
(monomers) and capsid intermediates can be related by 
the law of mass action ll34ll 



Pn<J^ = (picr^)"exp(-/3AG'„) 



(12) 



where p„ is the number density of an intermediate with n sub- 
units, a is the molecular dimension, /? = l/k^T is the in- 
verse of the thermal energy, and A(S'„ is the driving force to 
form an intermediate of size n. The driving force for assem- 
bly comes from the fact that subunits experience a favorable 
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FIG. 5: Examples of the influence of system parameters on assembly dynamics for design B3. (a) The fraction of complete capsids versus 
time, /c, is shown for 6m ~ 0.5 and Co ~ Na'^ /L^ — 0.11 at varying Eb illustrating the sigmoidal shape of capsid yields. Note that variations 
of /c are in discrete units of 0.06 because there are 1000 subunits and each complete capsid has 60 subunits. Variation of the final mass fraction 
of complete capsids, /c, is shown in (b)-(d): (b) variation with Co = 0.11 and 6^ ~ 0.5, (c) variation with Co at several values of Sb with 
6'm — 0.5, and (d) variation with 9m for several values of Sb and Co ~ 0.11. Note that eb does not denote the free energy to bind; there is a 
significant entropy penalty, calculated in Eg. 1151 



energy, £b, upon binding, but subunits also face an entropic 
penalty, which depends on the number of bonds and the local 
bonding network. 

The free energy for making a single bond to form a 
dimer can be determined by calculating the ratio of the par- 
tition functions for two bound subunits and two free subunits 



with 



92 /ql 



1 
8^ 



(13) 



exp [-(3u{Ri, R2, {bi}, {b2})] Hil, 2) 



where u is defined in Eq.|2l Jl2 describes the Euler angles of 
subunit 2, which specify the set of bond vector orientations, 
{b2}, and H{1, 2) is unity when u(Ri, R2, {bi}, {b2}) < 
—2ksT, and zero otherwise. In other words, we define two 
capsomers as bound if their potential energy of interaction is 
lower than —2ksT. The free subunits are taken to be at a 
standard state with unit density and free rotation, and the co- 
ordinate system is centered on Ri. 

Expansion of 2) to quadratic order in each coordinate 
about the minimum in the potential gives 



AG2 = -kBTliiq2/ql = -Eb - Tsb 



(14) 



3, 27:f3dW,{r) 
-2^" 



r=0 



1 2Pey 

2 eicbi 



(15) 



where the two terms represent translational and rotational en- 
tropy, respectively. This result is compared to binding free 
energies calculated with Monte Carlo simulations in Fig.|3] 

While there are many possible capsid structures consistent 
with most larger values of n, there is only one structure con- 
sistent with a complete capsid, which has n = N^^ subunits 
(Nc = 60 for the capsids studied in this work). The fact 
that misformed capsids and intermediates are generally not 
observed implies that AG„ is sharply peaked at n = iVc; de- 
fects that lead to larger or smaller capsids are unfavorable. 
There is a threshold density, pcc, at which the fraction of sub- 
units in capsids becomes significant IbI IstIIs^Is^ 

\np,,a^ ^ (3AGnJN,. (16) 

By analogy with Eq.O the free energy of a complete capsid 
can be written as 



AG 



-iVcnbeb/2-T(7V,-l)sbK), (17) 



where s\,{n\,) is the entropy penalty for a subunit in a complete 
capsid, where each subunit has n\, bonds. If we neglect the 
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dependence of the entropy penalty on the number of bonds 
[i.e. assume Sb("-b) ~ we can use Eqs.^][^ and 1 171 
to calculate pcc- Values of pcc for capsid design B3 {n\, = 3) 
with (pi-n = TT (the value used for all dynamical simulations 
with this work) are shown in Fig. |3 and are compared to 
kinetic assembly results in Fig.|6l 



IV. KINETICS OF CAPSID ASSEMBLY 

Capsid formation rate curves are sigmoidal. We have 
considered capsid assembly dynamics for design B3 (see 
Figs. 1 and 2) over ranges of subunit concentrations, Co (re- 
ported in dimensionless units, Cq = Na^/L^), binding en- 
ergies, Eb, and maximum binding angles 6^. The results we 
present use — tt; the effect of varying is similar to, but 
less dramatic than that of varying 0,,,. Dynamics of different 
capsid designs are discussed below. 

The fraction of subunits in completed capsids, /c, is shown 
as a function of time for several binding energies in Fig.|5t. 
In all cases for which significant assembly occured, the rate 
of capsid formation has a roughly sigmoidal shape. This is a 
general feature of assembly reactions \2($ that can be under- 
stood as follows. There is an initial lag phase during which 
capsid intermediates form and progress through the assembly 
cascade, followed by a rapid growth phase during which these 
intermediates assemble into complete capsids. Finally, growth 
slows when monomers (free subunits) are depleted and the 
remaining capsid intermediates are unable to bind with each 
other. 

Final capsid yields are non-monotonic with respect to 
parameter values, but high yields are possible. The frac- 
tion of subunits in complete capsids, /c, at the final observa- 
tion time, t — 6 X 10^, is shown in Fig |5] As Cq, Sb, or 
6'm increase, intermediates form and grow more rapidly, and 
thus capsid yields increase to as high as 90%, meaning 15 of 
the 16 possible capsids were completed. One of the primary 
results of this study is that a particle model that does not in- 
clude details such as heterogeneous nucleation or conforma- 
tional changes can predict such high capsid yields. 

Although capsids form more quickly as parameter values 
are increased, saturation of growth also occurs sooner and 
capsid yields are non-montonic in each parameter. The sensi- 
tivity of capsid yields to parameters seen in Fig. |5]is further 
illustrated with a kinetic phase diagram in Fig. |6] It demon- 
strates the coupled dependencies of capsid yields on system 
parameters. Phases are partitioned according to whether or not 
there is significant assembly, arbitrarily chosen as > 30%. 

The non-monotronic variation of capsid yields with pa- 
rameter values arises due to competition between faster 
capsid growth and kinetic traps. The initial steps in the as- 
sembly cascade result in the formation of fewer bonds than 
later steps. If the attractive energy of these bonds is not suffi- 
cient to overcome entropic loss the initial steps are uphill on 
a free energy barrier and hence, are slow. For parameter sets 
that are near p^c (see Fig. 13, where half of the subunits are 
in complete capsids at equilibrium, the fact that a complete 
capsid has many more bonds than initial assembly products 



indicates that the free energy barrier must be many times the 
thermal energy, k^T. Significant assembly, therefore, does 
not occur within the finite assembly times we consider un- 
til parameter values are much higher than the thermodynamic 
critical values, and we identify a kinetic lower critical surface 
(LCS) in Fig.|6|that bounds the regions with significant assem- 
bly from below and to the left. We consider results at finite 
observation times because capsid assembly reactions are lim- 
ited in-vivo by proteolysis times and in-vitro by experimental 
observation times. 

Increasing parameters increases the overall rate of capsid 
growth: higher subunit concentrations, Co, result in more 
frequent subunit collisions, higher values of 0^ increase the 
likelihood of binding upon a collision, and higher values of 
the binding energy, e^, decrease the rate of the reverse reac- 
tion (subunit unbinding). As parameter values cross the LCS, 
faster capsid growth leads to significant capsid yields, as seen 
m Fig.m At even higher values, however, assembly becomes 
frustrated by two kinetic traps (see Fig. and we identify an 
upper critical surface (UCS) in Fig|6lto the top and right of the 
regions in which assembly is kinetically accessable. Because 
of these kinetic traps, assembly only occurs at subunit-subunit 
binding energies that are much smaller than values calculated 
from atomistic potentials in Ref. 21 (see Table 1 of that Ref.). 
When the binding entropy (see Eq. I15> is included, however, 
the resulting free energies are consistent with association con- 
stants fit to assembly experiments with Hepatitis B capsids in 
Ref. 

We present results at three observation times to show how 
the distance between assembly boundaries expands in all di- 
rections as time increases. The rough boundary of the kineti- 
cally accessable region in Fig.|6lis a measure of the statistical 
uncertainty that results because each data point describes a 
single stochastic trajectory. For trajectories run with differ- 
ent random number seeds at a given set of parameter values, 
the final number of complete capsids typically did not vary 
by more than one capsid; however, variations during the rapid 
growth phase were larger. 

The kinetic trap depicted in Fig arises when progress 
through initial assembly steps is too rapid, allowing so many 
capsids to initiate that the pool of monomers and small in- 
termediates becomes depleted before a significant number of 
capsids are completed. If the remaining partial capsids have 
non-complementary geometries, further binding can only pro- 
ceed upon disassembly. This kinetic trap has been seen in ex- 
periments 1 12] and predicted theoretically |17]; this theory, 
however, assumes that only monomers can add to partial cap- 
sids. As discussed below, binding of capsid intermediates is 
an important mode of assembly and growth does not become 
frustrated until only intermediates with non-complementary 
geometries remain. 

The kinetic trap just described may not limit assembly in 
vivo, where there is a continual supply of new capsid pro- 
teins. Subunit bonding in configurations not consistent with a 
complete capsid (misbonding), however, can lead to a kinetic 
trap that could frustrate assembly even with an unlimited sup- 
ply of subunits. It is not surprising that misbonding occurs 
more frequently as 9^ increases, since there is a smaller driv- 
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FIG. 6: Changing model parameters reveals the kinetic phase diagram for design B3. Filled points denote parameter values for which 30% 
of subunits are in complete capsids (/c > 30%) by the observation time, tabs, open points indicate parameter values for which /c < 30%, 
and the dashed lines indicate the location of the thermodynamic critical surface, calculated with Eg. 1161 The first, second, and third columns 
correspond to observation time t„bf. ~ 1-2 x 10*, 8 x 10*, and 6 x 10''', respectively. The top row (a) shows cross-sections through Co and Sb 
with 9m = 0.5. The second and third row show cross-sections through and eb, with (b) Co — 0.11 and (c) Co ~ 0.037. In each case the 
region covered by filled points roughly defines a cross-section of parameter space within which assembly is kinetically possible. Simulation 
snapshots corresponding to the o in the right hand panel in row (b) and the A in the right hand panel in row (c) are shown in Fig.|7| 



ing force toward the minimum energy orientation. As noted 
by Berger and co-workers i23l . though, increasing concentra- 
tion and binding energy can stabiUze subunits with strained 
bonds, and the higher rate of capsid growth under these con- 
ditions can cause misbonded subunits to become trapped in a 
growing capsid by further addition of subunits. Because so 
many assembly pathways that do not lead to complete capsids 
are available at these parameter values, the minimum energy 
configuration, with complete capsids, is seldom realized and 
misformed capsids with spiraling or multi-shelled configura- 
tions dominate, as shown in Fig. Eb- Progression from this 
state to a completed capsid is extremely slow because break- 
age of many bonds is required. It would be difficult to assess 
the importance of configurations such as that shown in Fig. 05 
with models that have assumed assembly pathways. 

B5 capsids grow primarUy through additions of individ- 



ual subunits, while combination of clusters is essential for 
assembly of B4 and B3 capsids. The variations of final capsid 
yields with Sb for the capsid designs B3, B4, and B5 (see Fig 
13 are shown in FiglS^. Although assembly occurs within dif- 
ferent ranges of Sb for each design, assembly kinetics within 
these ranges are similar, as shown in FiglSji. In addition, the 
optimal assembly (/c « 0.9) for each design occurs at ap- 
proximately the same value of n\,eb ~ 50, meaning that the 
complete capsids all have about the same stability. Variation 
of assembly with Cq and 0^ (not shown) is also similar for 
each design. 

Although the capacity to assemble spontaneously is similar 
for each capsid design, the mechanism of assembly (near opti- 
mal assembly) for B5 is qualitatively different from that for B4 
and B3. We evaluated assembly mechanisms from simulations 
by tabulating the size, ribe, of the smallest intermediate in- 
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FIG. 7: Snapshots corresponding to unassembled points in Fig. |6| which illustrate the two kinetic traps described in the text, (a) Parameter 
values corresponding to the o in Fig. |6|allowed rapid capsid initiation and growth, which depleted monomers and small intermediates before 
capsids were completed, (b) At parameter values corresponding to the A, strained bonds are incorporated into growing capsids. These 
snapshots have zoomed in on a representative portion of the system to aid visibility. Subunit colors indicate the number of bonds: white = 0, 
red = 1, turquoise = 2, and dark blue = 3 bonds. 



volved in each binding or unbinding event (see Appendix B). 
The net contribution to capsid growth by intermediates of size 
nbe, &net("-be), IS shown for cach capsid design in Fig|9] While 
about 33% of subunits assembled as multimers (ribe > 1) for 
B3 and B4 capsids, multimer binding accounted for only 6% 
of all binding for B5 capsids. 

The influence of capsid design on mechanism can be under- 
stood by examining assembly pathways that are available if 
growth occurs only through monomer additions, such as those 
shown in Fig.[lO] Once dimerization occurs for B5, all sub- 
sequent monomers can add in such a way that two or more 
bonds are formed. For parameters values at which optimal as- 
sembly occurs, the formation of a single bond is unfavorable 
due to entropy loss, but the formation of two or more bonds is 
favorable. Consequently, the approximate projection of free 
energy onto cluster size (see Appendix C) shown in Fig. Illl is 
monotonically decreasing after formation of a dimer. 

For architectures B3 and B4, it is not possible to construct 
an assembly path for which monomers form multiple bonds 
at all cluster sizes greater than two (see Fig. I10> . Therefore, 
free energy profiles consistent with these architectures, shown 
in Fig.^] have numerous free energy barriers and local min- 
ima. At parameter values that are optimal for B5 assembly, 
progression of B4 and B3 intermediates through the assembly 
cascade is stymied by these free energy barriers. 

All free energy barriers vanish if parameters are increased 
until formation of single bonds is favorable. Since dimers are 
stable under these conditions, however, too many capsids ini- 
tiate and the system becomes mired in the kinetic trap shown 
in Fig. At moderate parameter values, B4 intermediates 
that correspond to local minima, such as trimers and hexam- 



ers, are metastable. Optimal assembly occurs when these in- 
termediates bind to growing capsids as fast as they are formed. 
The first local minimum for B3 capsids does not occur until a 
cluster size of five; in this case, optimal assembly occurs at 
parameter values for which dimerization is only slightly unfa- 
vorable and binding events involving dimers are common (see 
Fig.|9j. As parameter values are increased beyond optimal as- 
sembly conditions, formation of small intermediates becomes 
more rapid and multimer binding becomes more important for 
all capsid designs. 

Note that these free energies compare the relative stability 
of different multimers. There are also free energy barriers not 
shown that are associated with subunit binding or unbinding, 
which is required to transition between these states. 

Given the importance of multimer binding to assembly of 
some capsid designs, it might seem surprising that a rate equa- 
tion model ITtI Eoll that assumes capsids grow only via ad- 
dition of monomers predicts significant assembly for some 
parameter values. That particular class of models assumes, 
however, that monomers bind to intermediates with an aver- 
age binding free energy, independent of the number of inter- 
subunit contacts formed in a particular configuration. While 
not explicitly described herein, we have considered a model 
in which binding free energies were specifically calculated 
for each capsid intermediate, but the assumption that only 
monomers could bind to these intermediates was retained. 
The assembly dynamics predicted by this model for B5 were 
consistent with those seen in Brownian dynamics simulations, 
but concentrations of intermediates at local free energy min- 
ima built up for B4 and B3 designs, as well as for the de- 
sign shown in Fig IB of Ref . .20^ Because these intermediates 
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FIG. 9: B5 capsids grow primarily througli additions of individ- 
ual subunits, while binding of multimers to growing capsids is es- 
sential for assembly of B4 and B3 capsids. The fraction of bind- 
ing, 6nct (see Appendix B), for each cluster sizes ribc is shown 
for each design, and the total binding contribution for multimers, 
&muit = X]ni,.>i ^nairibt), is shown in the inset. Statistics were mea- 
sured for the parameters used in Fig. |8|through t = 2.25 x 10^, by 
which time the majority of assembly was completed for these param- 
eter values. 



FIG. 8: Assembly kinetics are similar for each capsid design, B3, B4, 
and B5. (a) The final capsid yields, /c at tobs = 6 x 10^, are shown 
versus binding energy for each capsid design for Co = 0.11 and 
6m = 0.5. (b) The assembly time series are shown for the optimal 
values of eb from (a) (Sb = 16.0 for B3, 12.7 for B4, and 10.5 for 
B5). 



could not bind with each other, the predicted assembly was 
much less efficient than that predicted by Brownian dynamics 
simulations. These results suggest that it may be important 
to consider binding of complexes that are larger than the ba- 
sic assembly unit when hypothesizing assembly mechanisms 
from experimental data. 

For many viruses, the basic assembly units are believed to 
be small intermediates, such as dimers or trimers. Experi- 
mentally observed high concentrations of these species dur- 
ing assembly reactions imply that they form rapidly and are 
extremely metastable. This feature could be included in our 
model by designing a new "subunit" that represents the basic 
assembly unit, or by choosing higher binding energies for cer- 
tain bond vectors. In this work, where all maximum binding 
energies are equal, the importance of multimer binding for B4 
and B3 capsids is not a result of the interaction between in- 
dividual subunits. Rather, the collective interactions of many 
subunits lead to a free energy profile with numerous local min- 
ima, which forces assembly to proceed through binding of 
intermediates. Although binding of trimers and other multi- 
mers is essential to B4 assembly in our simulations, the frac- 
tion of subunits that comprise these intermediates is always 
small compared to the fraction of subunits that are monomers 
or in completed capsids (see Fig. I12> . Hence, the signifi- 



(b) 

FIG. 10: (a) Examples of the assembly pathways that are available 
if only monomers can add to growing capsids. Once a B5 dimer is 
formed, subsequent monomers can always form two or more bonds, 
whereas all assembly paths for B3 and B4 require formation of single 
bonds, (b) One example of a multimer binding step for design B4 that 
avoids formation of single bonds. 



cance of binding of multimers would be difficult to detect by 
bulk experiments, such as light scattering or size exclusion 
chromatography, alone. The combination of these techniques, 
though, with selective deletion of residues through mutation 
|41| or single molecule experiments may offer insights into 
the importance of various intermediates in capsid assembly. 

Capsids are metastable in infinite dilution and capsid 
disassembly shows hysteresis. Since even single virions can 
sometimes infect cells, viral capsids must be metastable in in- 
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FIG. 11: Free energy projections onto the number of particles in 
a capsid, n, for each design, relative to monomers in solution at 
Co = 0.11 with 9m ~ 0.5. The values of Sb are the optimal val- 
ues from Fig.|8| 16.0, 10.5, and 12.7 for the B3, B4, and B5 designs, 
respectively. Free energies were obtained with Boltzmann weighted 
sums over configurations with a given cluster size as described in 
Appendix C. Lines are drawn between the points as a guide to the 
eye. 
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FIG. 12: The mass fraction of monomers, trimers, and complete 
capsids (60 subunits, with nt, — 4 bonds each) as a function of time, 
t, for design B4, illustrating that intermediate concentrations are al- 
ways small during successful assembly. The parameters are those 
given in Fig|8l 



finite dilution. Model capdsids also displayed this feature; for 
example, significant assembly with Cq = 0.008 did not occur 
for £b < 20.0 (see Fig.|6|l, but isolated complete capsids simu- 
lated without periodic boundary conditions (to model infinite 
dilution) were metastable through t = 6 x 10^ for Sb > 14.0. 
This result indicates that the final yield of the capsids at some 
binding energies will be different for a trajectory that started 
with mostly complete capsids then for a trajectory that started 
with random subunit configurations, as shown in Fig ^] In 
other words, there would be hysteresis between assembly and 
disassociation. 

Hysteresis arises because, as discussed above, there can be 
large free energy baiTiers separating the early stages of as- 
sembly, where there are few bonds per subunit, and complete 
capsids, for which each subunit has n\, bonds. Therefore, 
monomers (or complete capsids) can be metastable through- 
out a finite length trajectory above (or below) p^c- Once the 
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FIG. 13: Capsid yeilds as a function of eb (or equivalently inverse 
temperature or inverse denaturant concentration) illustrating hystere- 
sis between association and dissociation of capsids. The final cap- 
sid yields, /c at fobs = 6 x 10*, are shown for simulations started 
with subunits in random configurations (association) and simulations 
started with the final configuration from the B3 simulation shown in 
Fig.|8| which had /c — 0.9 (dissociation). Parameter values were 
Co = 0.11 and6i„ = 0.5. 

first subunit is removed from a metastable complete capsid, 
neighboring subunits have fewer bonds and further disassem- 
bly is rapid. Hysteresis in capsid assembly-disassembly has 
been seen with experiments and theory by Zlotnick and co- 
workers iTsli . 

V. CONCLUSIONS AND OUTLOOK 

In this work we examined ensembles of assembly trajec- 
tories for models that assemble into capsid-like objects. We 
demonstrated that computational models can distinguish driv- 
ing forces and corresponding mechanisms that lead to suc- 
cessful assembly from those that engender dynamic frustra- 
tion. Therefore, the approach we have outlined might be used 
to good effect to analyze the dynamics of other assembly mod- 
els (e.g., 1^ i23j 2^ 2M)- as well as that for other types of 
assembly models. In addition to the calculations presented 
herein, there are other ways that ensembles of assembly tra- 
jectories can be carefully analyzed. For instance, distributions 
of capsid formation times could be studied by simulation and 
potentially estimated with single molecule experiments; com- 
parisons could shed light on assembly mechanisms, such as 
the mechanisms for B4 and B5 assembly illustrated in this 
work. Trajectories generated by the approach described in this 
work can be a starting point for p erforming importance sam- 
pling in trajectory space |42, 43, 44], which would facilitate 
statistical analysis of ensembles of assembly pathways. 

Rechtsman and co-workers describe an "inverse statistical 
mechanical-methodology" that allows importance sampling 
in model space to design potentials that direct assembly 
into a particular ground state |45]. The ability to generate 
ensembles of dynamical trajectories invites a related strategy, 
in which one seeks to optimize a function of entire assembly 
paths, such as capsid assembly times or identities of key 
intermediates. This approach would involve importance 
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sampling steps in trajectory space, such as shooting moves 
146,1 . as well as sampling steps in model space. Understanding 
model features that lead to specific assembly behavior could 
guide and interpret experiments that involve mutated capsid 
proteins. 
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APPENDIX A 

The equations of motion given in Eq. |S]were integrated 
as follows. Translational displacements are calculated as de- 
scribed in Eq. 7 of Ref.js^ Bond vector orientations are spec- 
ified in body fixed coordinates at the beginning of the simula- 
tions. The space fixed coordinates, {b'^"'(i)}, are determined 
from a rotation matrix, A{t), which is evolved in time using 
quaternions, which satisfy the equation of motion given in Eq. 
3.37 of Ref. 31. This equation requires angular velocities, u), 
which are determined in an analogous fashion to the transla- 
tional displacements 

io, = 7,-/2 « + T^) + Sn (Al) 

where the torques are calculated at two points 

Tf(t) EE r,({R,(t),b(t)}) 

r^{t) ^ T,{{^\{t),h\m (A2) 

with the predictor positions determined as in Ref. 

R,' = R, + (5t(7F, + ^F,). (A3) 

The predictor bond orientations, {b''(t)}, are determined from 
a predictor rotation matrix, which is calculated from Eq. 3.37 
of Ref. 3 1 , using predictor angular velocities calculated as 

i^\^l.Tt + 5n. (A4) 

This formulation assumes that subunits are hydrodynamically 
isolated and that rotational and translational sources of friction 
are not coupled; these assumptions can be relaxed as in Ref. 

APPENDIX B 

The contribution of multimer-binding to the final assembly 
product was calculated from simulations as follows. Multi- 



mers were designated by clustering subunits connected by one 
or more bonds. A binding event occurred when the size of a 
cluster changed, either through a combination of two clusters 
(positive binding) or division of a cluster (negative binding). 
Cluster sizes were output every 10 steps; more than one bind- 
ing event involving the same cluster within 10 steps was found 
to be exceptionally rare. The size of a binding event, ribe, was 
defined by the size of the smallest reactant cluster for positive 
binding or the smallest product cluster for negative binding. 
The net forward binding due to events of size 7ibe is given by 

6net(?lbe) = ^^be (fc+(»^be) - 6_(nbe)) (Bl) 

where 6+ and &_ are the number of positive and negative bind- 
ing events of size ribe, respectively. 



APPENDIX C 

The free energy, Ag^, to build a particular capsid configu- 
ration, i, from a bath of subunits at concentration Co, can be 
determined by analogy to Eq.[2]if the dependence of binding 
entropy on the number of bonds is neglected 



^9. = J2 [-&jeb/2]-(n,-l)T [fcB ln(7ra3Co/6) + .Sb(l)] , 

(CI) 

where j sums over each subunit in the capsid, rii is the number 
of subunits in configuration i, and bj is the number of bonds 
for subunit j. We project the free energy onto the number 
of particles in a capsid, n, by summing over configurations 
containing n of subunits 

G„ = -fcBTln V ,5„.,„ exp (-/3AgO . (C2) 

We performed this sum with Monte Carlo simulations in 
which trial capsid configurations were generated by adding 
or deleting subunits from current configurations, and then ac- 
cepted or rejected according to the Metropolis |48| criterion, 
with the Boltzmann distribution given by Eq. lCll Only config- 
urations consistent with minimum energy bonding were con- 
sidered; thus, this approach was only applicable for parame- 
ters with which misformed capsids do not occur. The average 
free energy for capsids of size n was efficiently calculated by 
carrying out umbrella sampling |34] in which a harmonic po- 
tential as a function of capsid size was used to bias the number 
of subunits in the capsid. 
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